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Abstract 

Among the 13 TLRs in the vertebrate systems, only TLR4 utilizes both Myeloid differentiation factor 88 (MyD88) and Toll/ 
lnterleukin-1 receptor (TIR)-domain-containing adapter interferon-fi-inducing Factor (TRIF) adaptors to transduce signals 
triggering host-protective immune responses. Earlier studies on the pathway combined various experimental data in the 
form of one comprehensive map of TLR signaling. But in the absence of adequate kinetic parameters quantitative 
mathematical models that reveal emerging systems level properties and dynamic inter-regulation among the kinases/ 
phosphatases of the TLR4 network are not yet available. So, here we used reaction stoichiometry-based and parameter 
independent logical modeling formalism to build the TLR4 signaling network model that captured the feedback regulations, 
interdependencies between signaling kinases and phosphatases and the outcome of simulated infections. The analyses of 
the TLR4 signaling network revealed 360 feedback loops, 157 negative and 203 positive; of which, 334 loops had the 
phosphatase PP1 as an essential component. The network elements' interdependency (positive or negative dependencies) 
in perturbation conditions such as the phosphatase knockout conditions revealed interdependencies between the dual- 
specific phosphatases MKP-1 and MKP-3 and the kinases in MAPK modules and the role of PP2A in the auto-regulation of 
Calmodulin kinase-ll. Our simulations under the specific kinase or phosphatase gene-deficiency or inhibition conditions 
corroborated with several previously reported experimental data. The simulations to mimic Yersinia pestis and £ coli 
infections identified the key perturbation in the network and potential drug targets. Thus, our analyses of TLR4 signaling 
highlights the role of phosphatases as key regulatory factors in determining the global interdependencies among the 
network elements; uncovers novel signaling connections; identifies potential drug targets for infections. 
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Introduction 

Toll-Like Receptors (TLRs), thirteen in vertebrates, are the 
members of pattern recognition receptor family (PRRs) and 
recognize the pathogen-associated molecular patterns (PAMPs) 
[1,2]. Upon ligand binding TLR signal is processed in a Myd88- 
dependent and MyD88-independent (TRIF dependent) manner 
[3], where MyD88 and TRIF are the adaptor molecules 
differentially recruited to the TLRs. Among all the TLRs, only 
TLR4 utilizes both Myd88-dependent and TRIF-dependent 
pathways [3]. This makes TLR4 signal processing a relatively 
complex process as compared to other TLR family receptors. 

In such complex biological systems of signal processing and 
integration, understanding the inter-regulation between the 
signaling elements by computational reconstruction of the 
signaling networks enables one to systematically investigate the 
regulatory principles associated with the network and thus to build 
hypothesis that can subsequently be tested through experiments. 
For capturing the quantitative system level properties of signaling 
network dynamic modeling approaches [4,5] are best suited 
provided kinetic details of the interactions and the concentrations 
of the species of the network are known. But for signaling networks 



as large as of TLR4 knowledge of kinetic parameters/ concentra- 
tions of the elements in the model is extremely sparse; nonetheless, 
the interaction details among the elements of the TLR4 network 
are relatively well characterized. Such knowledge of reaction 
stoichiometry was used for building qualitative Boolean logic 
based models independent of kinetic parameters or concentration 
of the network components [6]. By adopting such logical modeling 
formalism, which utilizes the available molecular interaction 
details and reaction stoichiometry of the network to convert the 
signaling interactions to logical connections [6,7], we have 
analyzed here the TLR4 signaling networks. 

Here, we have constructed a logical model of TLR4 signaling 
utilizing information derived from hundreds of independent 
experimental reports. The preliminary information was gathered 
from the comprehensive map of TLR signaling [8] which we 
further updated with recent experimental findings. The experi- 
mental information was converted to logical connections [7] . The 
model identified a total of 360 operational feedback loops of which 
most of the loops (positive or negative) owe their origin to four 
phosphatases, MKP-1, MKP-3, PP1 and PP2A. We calculated the 
relative contribution of each phosphatase in determining the total 
number of feedback loops in the system by systematically carrying 
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out several in-silico knockout studies and subsequently analyzing 
the changes in global dependency (positive/negative/neutral) 
among the network elements. The model was useful in extracting 
previously unreported/less emphasized paths of signal propagation 
for providing plausible explanations to paradoxical experimental 
results. For example, the counteractive roles for ERK-1/2 in the 
augmentation or inhibition of IL-10 could be explained by 
extracting novel signaling pathways from the model. We tested the 
predictive power of the model by successfully simulating effects of 
various experimentally reported knockout conditions. Finally, we 
subjected the model to in silico perturbations that mimic the effect 
of infections by Escherichia coli and Yersinia pestis. The model at its 
current state is open to experimental tests and thus further fine- 
tuning; as it was observed that training such logical models with 
high throughput experimental data remarkably increases its 
predictive power [9]. 

Results 

Constructing the logical model of TLR4 signaling 

We started the model building by extracting the information on 
the TLR4-specific interactions from the comprehensive map of 
TLR signaling [8]; newer information (published later, after [8]) 
was subsequendy incorporated in the model (Table SI). All the 
reactions were next converted to logical connections (hyper- 
graphs), creating 181 species (nodes) and 263 reactions (edges) 
based on Boolean logic functions AND, OR, NOT and their 
combinations [6,10] (Figure 1). The logical functions were 
implemented in accordance to the nature of biological connections 
between the network elements using CellNetAnalyzer [10]. Table 
S2 and Table S3 show the logical reactions and species, 
respectively, involved in the TLR4 signaling. The nodes in 
signaling network are kinases, phosphatases, adapter molecules, 
transcription factors, cytokines etc. The edges are the unidirec- 
tional hypergraphs, which represent the direction of signal flow. 

In the logical model the state 1 and 0 represents on (active) and 
off (inactive) state, respectively. 

In the simplest kind of interaction where an upstream directiy 
activates its downstream, activation state 1 of the upstream 
activator is transferred to the downstream, representing the signal 
flow. Consider for example, a set of molecules X, Y and Z; 1] if 
activation of Z requires activation of both X and Y simultaneously, 
then the logic gate AND connects A, B to C as: X AND Y = Z; 2] 
If Z is activated when any of its upstream X or Y are active, then 
the Boolean gate OR connects X, Y, Z as: X OR Y = Z; 3] If Z is 
active only when X is active but Y is NOT active then the NOT 
gate connects the elements as X+NOT Y = Z. Figure 2 (methods 
section) shows a set of representative molecules from the 
constructed TLR4 network whose biological interactions required 
representation utilizing AND, OR and NOT gates. 

The TLR4 model comprised of 40 input nodes. These inputs 
correspond to TLR4 receptor; TLR4 ligand, co-receptor CD 14 
and lipopolysaccharide binding protein (LBP) and the rest 36 are 
the intermediates for which direct activators or inhibitors are not 
know clearly. Table S2 shows the list of all input nodes used in the 
model. The intermediate nodes carry the signal from the input 
layer which mostly comprises components of various evolutionarily 
conserved kinase-signaling pathways. Complexity of the signaling 
increases as the signal goes down from the membrane to the 
cytoplasm and then to the nucleus as more and more intermediate 
signaling components is activated by the incoming signal. Finally, 
the signal converges to output elements comprising various pro- 
and anti-inflammatory cytokines, transcription factors and various 



target genes, trigger of cell division, transcription factors, apoptosis 
and regulation of the RNA metabolism. 

In biological system activation of signaling elements in a 
pathway is temporal in nature which is determined by the 
hierarchy of signal flow from membrane to nucleus. To mimic the 
dynamicity (or pseudo dynamics) of signal flow like a real system 
we have opted time scale scenario formulation by selectively 
assigning different time scale to different signaling elements in a 
hierarchical manner [11]. We implemented time scales as below: 

Time scale 0: Formation of ligand receptor complex and 
activation of housekeeping species. 

Timescale 2: Activation of the all-cytoplasmic reactions taking 
place after recruitment of ligand to the receptor. 

Time scale 4: Activation of reactions corresponding to the nuclear 
translocation of specific cytoplasmic kinases that leads to triggering 
of transcription factors or activation of other target molecules/ 
genes. 

Time scale 7: In this time scale signaling events in the nucleus that 
are triggered by the MyD88 dependent pathways were activated. 
List of the reactions involved in the MyD88 dependent pathway is 
given in Table S2. 

Time scale 7.1: Activation of the TRIF dependent pathway was 
assumed to occur at later as suggested by experimental studies 
[12]. List of reaction involved in the TRIF dependent pathway is 
given in Table S2. 

Time scale 8: Finally in the time scale 8 the cytoplasmic and 
nuclear phosphatases and inhibitors of the TLR4 signaling were 
activated. This leads to complete shutdown of information flow in 
the system. 

In the in-vivo settings, the cytoplasmic/nuclear phosphatases 
can be activated earlier to the time required for the signal to reach 
the nucleus or they remain constitutively active causing an 
attenuation/suppression of the signal strength, thus imparting 
dynamically achieved finer regulations [4,5,1 1,13], so activation of 
phosphatases attenuates the signal flow in most cases but may not 
always abolish it. However, owing to the nature of our modeling 
approach (output is either 0 or 1), such finer modulations of the 
processed signal cannot be captured [7,14]. Nonetheless using the 
logical models analyses can be performed on various structural 
aspects of regulations' in the signaling network. Here we focused 
our analysis on understanding role of the phosphatases in 
regulating the global interdependencies among the elements of 
the TLR4 signaling network. 

TLR4 signaling network contains a large network of 
positive and negative feedback loops 

Feedback loops are pivotal for the functioning of almost all the 
signal transduction networks. Various dynamic analyses of 
signaling systems show that positive feedback loops lead to signal 
amplification, multi-stability and hysteresis or tolerance [15,16] 
whereas negative feedback loops trigger signal attenuation, 
adaptation, oscillations or delay in activation [17,18]. From the 
large signaling systems with information on reaction stoichiometry, 
logical analysis can identify the feedback loops and their nature 
(positive or negative). A feedback loop- whereby a signal from a 
signaling intermediate "I" flows through the intermediate "L" and 
comes back to influence the activity of "I"- is an interaction graph- 
based property of the model [6] . So, during the calculations, the 
number of times such conditions are satisfied yielded the number 
of feedback loops in the system. 

An important element in the feedback loops are the phosphatase 
of the network [14]. Our simulations identify 360 feedback loops, 
of which 157 were negative and 203 were positive feedback loops. 
Most of these feedback loops emerged from the kinase-phospha- 
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Figure 1. TLR4 Signal Transduction pathway map. Logical model of the TLR4 signal transduction network built in CellDesigner. The input 
species are shown by arrow pointing to the species while output species are shown by arrow going away from a species. Final outputs of model are 
represented as yellow colored boxes. Black dot in picture indicates AND gate. Black colors of the edges represent the activatory signals to the 
downstream node while the inhibitions are represented in red color. 
doi:10.1371/journal.pone.0092481.g001 



tase interactions in the network. We explored the role of 
phosphatases involved in TLR4 signaling and calculated their 
contribution in the total feedback loops of the network. Figure 3 
shows the schematics of relative contribution of different 
phosphatases, alone and in combinations in determining the total 
number of feedback loops in the system. To decipher the relative 
contribution of various phosphatases we systematically knocked 
out a target phosphatase and calculated the remaining feedback 
loops for each of the knockout conditions (Figure 3). The 
phosphatase PP1, as we found from the in-silico knock down, is 
the most contributory phosphatase in determining the total 
number of feedback loops (Figure 3) as ~93% (334 out of 360) 
feedback loops required PP1 as an essential component. 

When all the phosphatases of the system were knocked out, 25 
feedback loops remained. These 25 feedback loops owed their 
origin purely to kinase-mediated interactions of which 5 were 
negative and 20 were positive feedback loops. The interactions 
that gave rise to the 25 feedbacks were next extracted and shown; 
1 1 reactions were found to comprise the 5 negative feedback loops 
and 31 reactions comprised the 20 positive feedback loops, 
respectively (Table S4). Relative contribution of each of the 
reactions in the 25 feedback loops was also calculated and shown 
(Table S4). Notably, 4 out of 5 negative feedback loops comprised 
NF-kB as an essential component and Vavl was found to be 
essential element in more than 60% of the 20 positive feedback 
loops. 



Interdependency among the network elements: wild 
type and phosphatases knockout conditions 

As interdependencies among the signaling intermediates char- 
acterizes cell signaling [9,19], any local changes imposed in the 
system can possibly lead to global alteration of interdependencies 
among the remaining elements in the network in varying 
degrees[9,l 1]. Experimental investigations and hypothesis testing 
commonly involve targeted knockout studies [20-22] although 
systems level implications of such knockout studies are less 
understood. As our analyses suggest (Figure 3), knockout of 
different phosphatases do change the number of feedback loops 
differentially, indicating that the interdependencies among the 
network elements plausibly changes with each perturbation. 

Typically in large-scale networks a remote downstream kinase 
could be an activator/inhibitor of an upstream molecule and the 
process of activation/inhibition comprises multiple intermediates 
[17,18]. Here positive or negative influence of a network element 
on the rest of the network elements is calculated using the 
dependency matrices [6,10]. For example, a species X can 
influence any other species Y of a network in various ways: X 
could be a total activator (TA), total inhibitor (TI), one of the 
inhibitors (I) or one of the activators (A). Also species X can have 
no influence (N), or have both positive and negative influence (U, 
called as ambivalent factor) on species Y, all of which could be 
calculated using the dependency matrix [6]. 

The model of the TLR4 signaling network contained 181 
species, so the size of the dependency matrix is 181x181; making 
its visual comprehension and representation inconvenient. We 
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Figure 2. Conversion of reactions from the comprehensive map to the logical reactions. Set of representative interactions extracted from 
comprehensive map of TLR signaling and represented as interaction graph and interaction hypergraph. 
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26 



\ 


L 

1 


168 





(ppi)-/- 






(MKPl)-/ 




} 


r 










r 


28 




25 




26 



(PP2A+MKP1)-/- (PP1+PP2A)-/- 



25 



(PP1+MKP1)-/- 




25 



(PP1+PP2A+MKP3)-/- (PP1+PP2A+MKP1)-/- 



(Wild Type) 



360 



25 



J 

} 


f 




210 



31 



(MKP3)-/- 

-L> 



(PP2A)-/- 



118 




31 




26 



(MKP3+MKP1)-/- (MKP3+PP2A)-/- (MKP3+PP1)-/- 



26 



28 



(PP1+MKP3+MKP1)-/- (PP2A+MKP3+MKP1)-/- 



(PP1+PP2A+MKP3+MKP1)-/- 



Figure 3. Relative participation of phosphatases in feedback loops. Figures show number of remaining feedback loops in different 
phosphatases knockout conditions. 
doi:1 0.1 371 /journal.pone.0092481 .g003 



PLOS ONE I www.plosone.org 



4 



April 2014 I Volume 9 | Issue 4 | e92481 



Logic-Based TLR4 Signaling Analyses 



thus selected a subset of the total elements of the complete 
dependency matrix comprising molecules such as the MAP kinases 
(ERK, p38mapk, JNK), Akt, NIK, IKK and CAMKII which are 
representative of distinct signaling pathways activated during 
TLR4 signaling, together with the phosphatases MKP-1, MKP-3, 
PP1, PP2A. The chosen kinases are also pivotal biologically, as 
they carry out growth, differentiation, and proliferation and are 
involved in several diseased conditions [23-28]. 

Effects of various in-silico phosphatases' knockout and resultant 
alterations of the positive/ negative dependencies among species 
are detailed (Table S5). The complete dependency matrices for 
wild type and various knockout conditions are also shown 
(Supporting information). A total of 16 conditions were simulated: 
wild type, MKP-3" 7 ", PP2A" 7 ", MKP-1" 7 ", PP1" 7 ", (MKP-3+ 
MKP-1)" 7 ", (MKP- 3+PP2 A) ~ 7 ~ , (MKP-3+PP1)" 7 ", (MKP-1 + 
PP2A)" 7 ", (MKP-1+PP1)" 7 ", (PP1+PP2A)" 7 ", (MKP-3+MKP- 
1+PP2A)" 7 ", (MKP- 3+PP2 A+PP 1 ) ~ 7 ~ , (MKP- 1 +PP2 A+PP 1 ) ~ 7 
" and (MKP- 3+MKP- 1 +PP2 A+PP 1 ) ~ 7 ~ . Several interesting 
behaviors were observed in the knockout conditions all of which 
cannot be discussed due to the volume of the information. So, we 
state some interesting regulatory behaviors exhibited by the system 
in different knockout conditions that can be easily tested 
experimentally. 

> ERK- 1/2 acted as an ambivalent factor of MEKK3 and 
IKK in the wild type condition but in the MKP-3 knock out 
condition, ERK- 1/2 become an activator of MEKK3 and IKK; 
this suggests that the negative effect from ERK- 1/2 to MEKK3 
and IKK is mediated in an MKP-3 dependent manner. In 
addition, in the MKP-3 knockout condition ERK- 1/2 becomes an 
inhibitor of PP2A and p38MAPK whereas in the wild type 
condition ERK- 1/2 is ambivalent to both. ERK- 1/2 becomes an 
activator of CAMKII and inhibitor of Akt in MKP-3 knock out 
condition whereas in wild type condition ERK- 1/2 is ambivalent 
to both (Figure 4A-E). 

> When PP1 is knocked out MEKK3 becomes an activator of 
PP2A and Akt and inhibitor of itself, IKK and CAMKII 
respectively, whereas in the wild type condition its influence is 
ambivalent for all of these kinases (Figure 5C). An experimentally 
testable hypothesis from this the analysis would be: MEKK3 
inhibits PP2A and Akt but activates itself, IKK and CAMKII 
(Figure 5A-E) in the wild type setting which could be tested through 
perturbation studies like overexpression and/ or si-RNA mediated 
inhibition of MEKK3 and compare the effect of such perturbation 
on these kinases with their respective wild type counterparts. 

> In PP2A knockout conditions, CAMKII becomes its own 
inhibitor and its influence on AKT is completely abolished (N) 
whereas in the wild type condition CAMKII ambivalently affects 
itself and Akt. Thus the positive self-regulatory loop of CAMKII 
must have PP2A as an essential element. Also PP2A has a dual role 
(positive and negative) in CAMKII mediated regulation of Akt in 
the wild type condition. Further, knockout studies in the model 
show that the positive loop from CAMKII to Akt is abolished 
when MKP-1+PP1 is knocked out (Figure 6A-B). 

> All the signaling paths from MKP-1 to MEKK3, CAMKII 
and Akt have PP2A as an essential element as in the PP2A knock out 
condition MKP-1 is no longer connected to these kinases 
(Figure 7A-C). Such coherent regulation between the phosphatases 
can also be tested experimentally by systematic knock down studies. 

Coexisting signaling paths show plausible positive and 
negative pathways for ERK-1/2 mediated IL-10 
production 

Independent experiments show that in TLR4 signaling ERK- 1 / 
2 regulates IL-10 production both positively [29,30] and 



negatively [31]. As analysis of logical models show that complex 
connectivity among the elements of large networks lead to 
ambivalent effect [14], we explored whether both positive and 
negative signaling paths between ERK- 1 / 2 and IL- 1 0 potentially 
exists in the wild type condition. Due to the overwhelming 
complexity of connections in the network the total number of 
signaling paths from ERK-1/2 leading to IL-10 were found to be 
2508, out of which exactly 50% (1254) were positive and the rest 
50% were negative. Notably, all of these signaling paths contained 
PP2A and MKP-3 as essential components, because knockout of 
either of these phosphatases resulted in 0 signaling paths between 
ERK-1/2 and IL-10. Figure 8 shows the shortest paths that 
connect ERK-1/2 to IL-10 both positively and negatively. 
Similarly, we extracted a pathway that positively connects ERK- 
1/2 to IL-12, which is mediated through NF-KB (Figure 9) 
[32,33]. 

Logical steady state analysis and model validation using 
experimental results 

Logical steady state analysis (LSSA) was performed to under- 
stand the characteristics of signal flow through the TLR4 network. 
In LSSA the state 1 means an on state and a state 0 means off 
state. In LSSA, the state of each model element is consistent with 
the value of its associated Boolean function, and therefore, once a 
Boolean network has moved into a logical steady state, it will stop 
to switch and then retain this state [6] . The signal starting at input 
node (TLR4 receptor/ligand) traverse the entire span of the 
network through the series of logical gates and the cumulative 
result of such logical interactions decide the 0 or 1 state of each 
element of the network. Logical steady state was calculated at 
different time scales (biological logic behind the time scales 0, 2, 4, 
7, 7.1, 8 is explained in the model building sections and also given 
in supporting information). For example, we showed the LSS 
corresponding to time scales 0, 2, 4, 7, 7.1 in the wild type 
condition and subsequently effect of different knockout/knockin 
conditions at the time scale 7.1. At time point 7.1 (as explained in 
methods section) all the elements of the network are activated, 
hence global effect of perturbations was captured in time scale 7.1 
(Table S6). Simulation of the model with the next higher time scale 
(time scale 8) corresponds to complete shutdown of the signal in 
the system. 

As demonstrated earlier [14], results of LSS can be qualitatively 
compared to experimental results. Here we imposed various 
experimentally observed scenarios for TLR4 signaling and 
compared the simulated outcomes with the experimental reports. 
Table S7 lists several simulation results that corroborated with the 
experiments. When the simulation results matched the experi- 
mental results, we extracted the signaling path that links the 
perturbed species to the affected species. This was done because in 
experimental studies only the cause (perturbed species) and effect 
(changes in target species) are observed and the signaling path/ 
information processing route that connects both the species are 
often considered as a black box. For example, it was recently 
reported that in the TRADD 7 mice phosphorylation of the 
MAPKs' (p38MAPK, ERK-1/2, Jnk) and I-kB is reduced [34]; 
this is qualitatively observed in our model as 0 activation states of 
the MAPK's and 1-Kb when TRADD = 0. We found 6 positive 
signaling paths that connect TRADD to MAPKs' and 110 
signaling path that connects TRADD to I-kB. The signaling paths 
unravel the plausible intermediates involved in transmitting the 
signal from TRADD to the affected molecules, which is open to 
experimental validation. 
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Figure 4. Influence of ERK1 12 on various kinases in wild type and MKP3 knockout condition. A. ERK1 12 regulates MEKK3both positively as 
well as negatively in wild type condition. In the MKP3 knockout condition the regulation of MEKK3 by ERK1/2 is only positive; similarly the effect of 
the ERK1/2 on IKK (B), p38MAPK (C), Akt (D) and CAMKII (E), in wild type and MKP3 knockout condition is shown. Black arrow represents positive and 
red arrow represents negative influence respectively. 
doi:10.1371/journal.pone.0092481.g004 



Therapeutic utility of the model: Pathogen specific 
perturbations in the logical network and proposition of 
drug targets 

Several pathogens are known to perturb TLR4 signaling during 
their host invasion. We tested two cases of infection; infection by 
uropathogenic E.coli and Yersinia pestis and investigated whether the 
model predictions fall in the same lines as the experimental reports 
on the same. 

Yersinia pestis infection. Yersinia pestis is a Gram-negative 
bacterium that causes bubonic plague. YopJ protein found in 
Yersinia inhibits MKK6 and IKK by acetylation [35]. YopJ is also 
shown to inhibit the TRAF6 and TRAF3 by deubiqutinization 
[36]. Yersinia Yopj protein also inhibits MKK6 and IKK [37]. We 
computed the LSS for such pathogenic conditions by assigning 
zero value to the YopJ affected molecules (MKK6, IKK, TRAF6 
and TRAF3). LSS values of different species in the condition of 
infection are shown (File SI). LSS shows that pro-inflammatory 
cytokines production is inhibited while anti-inflammatory cytokine 
IL-10 remained uninhibited which was also observed experimen- 
tally [35,37]. Yersinia plausibly uses this strategy to overcome the 
defense mechanism of the immune system by selective inhibition of 
the pro-inflammatory cytokines' production and allowing only the 
production of anti-inflammatory cytokine such as IL-10. IL-10 
production in the infection condition is ERK1.2 dependent, which 
is activated through a path comprising Rac.GTPase. To alter the 



infection condition one strategy would be to suppress the ERK- 1 / 
2 activation which in our model would mean ERK-1/2 = 0. We 
thus calculated the MIS (Minimum Intervention Set [6], see 
methods section) for IL-10 production by assigning the MIS goal 
to be ERK1.2 = 0, while keeping MKK6, IKK, TRAF6 and 
TRAF3 = 0 such that the effect of Yersinia infection is considered as 
starting condition in the model. Table 1 shows a set of molecules 
that are calculated as MIS suitable for inhibiting IL- 1 0 production 
under Yersinia infection. 

Uropathogenic E.coli infection. Strains of Uropathogenic 
E.coli are responsible for the urinary tract infection [38]. They 
secrete Tcpc protein containing the TIR like domain that binds to 
Myd88 making it functionally ineffective [39]. In the model we 
have implemented the condition by assigning zero value to Myd88 
and then simulated the model to calculate the LSS (File SI). 
Experiments show that during Uropathogenic E.coli infection TNF-a 
production is abrogated [37,38]. Simulation of the model with 
Myd88 1 condition shows that TNF-a production is zero, even 
when the TRIF dependent pathways were considered activated. 
We performed MIS calculation by setting a MIS goal of 
reactivating TNF-oc in a Myd88 independent manner. We found 
363 MIS sets for this goal. Set of molecules whose coupled 
activation/inhibition would result in reactivation of TNF-a is 
given in Table S8. 
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Figure 5. Auto regulation of MEKK3 and its influence on different kinases/phosphatase in wild type and knockout conditions. 

Influence of MEKK3 on PP2A (A), Akt (B), MEKK3 (C), CAMKII (D) and IKK (E), in wild type and PP1 knockout condition is shown. Black arrow shows 
positive and red arrow represents the inhibitory influences respectively. 
doi:1 0.1 371 /journal.pone.0092481 .g005 



Discussion 

Here, we present the logical model of TLR4 signaling, 
comprising 181 nodes and 263 edges, which is to the best of our 
knowledge the largest logical model of signaling network built till 
date. Despite having incredibly complex interactions among the 
molecules of the network it was possible to extract several novel 
regulatory properties that owe their origin to stoichiometry of 
interactions and their deduced Boolean relations based on such 
stoichiometry matrix. Implementation of logical modeling meth- 
odology for large scale signaling networks is a relatively recent 
advancement in systems biology [6] . Till date, the prominent logic 



based models of large scale signal transduction are- T cell signaling 
network with 94 nodes and 123 interactions [7], Apoptosis 
signaling network with 86 nodes and 125 interactions [14], 
EGFR/ErbB signaling network with 104 nodes and 204 interac- 
tions [9] . All these models showed that logic based modeling can 
be of immense help in extracting a systems regulatory design 
stemming from the topology of interaction among its components. 

TLR4 signaling network comprises at least 181 nodes and 263 
interactions, making it one of the most complex signal processing 
systems to be subjected to quantitative modeling. In an earlier 
attempt, semi-quantitative model of TLR4 signaling was con- 
structed where effect of controlled perturbation of input on the 
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Figure 6. Self-regulation of CAMKII. (A) CAMKII and Akt (B) in the wild type and PP2A knockout out condition are shown. Black arrow shows 
positive influence; red arrow shows the inhibitory effect and cross mark shows no effect. 
doi:1 0.1 371 /journal.pone.0092481. g006 
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Figure 7. Influence of MKP1 on MEKK3, Akt and CAMKII. (A) MEEK3, (B), Akt and (C) CAMKII, in the wild type and PP2A knockout out 
conditions are shown. Black arrow shows positive influence; red arrow shows the inhibitory effect and cross mark shows no effect. 
doi:1 0.1 371 /journal.pone.0092481 .g007 



system's output was modeled [40], capturing properties of the 
network like redistribution of signal flow between the MyD88 and 
TRAM-dependent pathways upon inhibition of either of them. 
However, some of the crucial regulatory aspects of the TLR4 
signaling network such as identifying the feedback loops within the 
network, multi-species interdependencies in coherently determin- 
ing the output of the system was not addressed till now. 

We built the logic based model of TLR4 signaling and largely 
focused our analysis on understanding the role of phosphatases in 
determining the global emergence of feedback loops and 
interdependencies among the network elements. It was previously 



ERK-1/2 
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Nf-kB 



TTP 



IL-T1 



Figure 8. Shortest positive and negative signaling path from 
ERK1/2 to IL-10. Positive and negative signaling path from the ERK1/2 
to IL-10 is shown. Red bar show inhibition and black shows activatory 
signal. 

doi:1 0.1 371 /journal.pone.0092481. g008 



shown that phosphatases play crucial roles in the emergence of 
feedback loops in signaling pathways [41—45]. The in-silico studies 
performed here using knockout of various phosphatases indicate 
that phosphatases orchestrate a complex network of feedback 
loops, implying interdependencies among the phosphatases 
themselves; indeed, in a smaller network of MAPK signaling in 
B cells, the phosphatases' interdependency in regulating kinases' 
phosphorylation has been already implicated [45]. Because the 
phosphatases work in unison in the wild type system, it is hard to 
decipher the contribution of an individual phosphatase in 
influencing the interdependency among kinases in the network. 
We found that systematic knock out of one or multiple 
phosphatases in combinations differentially altered the interde- 
pendency between phosphatases, between kinases and between 
kinases and phosphatases. We systematically compared the effect 
of various phosphatase knockout conditions and studied the 
change in regulations among the pair of network elements whose 
positive or negative interaction is known for the wild type 
conditions. We show that the nature of regulations among various 
elements of the network differentially changes as a result of a 
phosphatase knockout. Although dynamic models and model- 
based experimental investigations [5,15,16,45,46] revealed the 
regulatory mechanisms underlying various experimentally ob- 
served phenomena in signaling pathways, global reorganization of 
the interdependency between network elements in a specific knock 
out condition remains poorly understood in general. Our analysis 
of the in silico knockout conditions presented several experimen- 
tally testable hypotheses that aim to elucidate the structure based 
design principles in the TLR4 signal in a finer way. 

We also showed the advantage of logical model formalism in 
extracting the hitherto unknown signaling pathways between the 
network elements. For example, the experimentally observed 
counter-regulation of TLR4-induced ERK-l/2-mediated IL-10 
production [29-31] is now shown to be due to the positive and 
negative signaling paths from ERK- 1 / 2 and IL- 1 0, dominance of 
either type of path could thus result from cell type specific or 
context specific expression/ activation of intermediates connecting 
ERK- 1 / 2 and IL- 10. However the nature of information flow here 
restricts one from inferring the relative dominance of either of 



PLOS ONE | www.plosone.org 



8 



April 2014 | Volume 9 | Issue 4 | e92481 



Logic-Based TLR4 Signaling Analyses 



ERK1/2 



MSK1/2 



NF-kB 



i 



IL12 



Figure 9. Shortest positive path form ERK1/2 to IL-12. Positive signal path from the ERK1/2 to the IL-12 is shown. 
doi:1 0.1 371 /journal.pone.0092481 .g009 



positive or negative connections between ERK-1/2 and IL-10. 
This is again because of the definition of on state (=1) and off state 
( = 0) of a species, as once ERK- 1/2 attains an on state it will 
transmit the information to IL-10 through intermediates irrespec- 
tive of wiring of signaling intermediates upstream to ERK-1/2. 
For example when we simulate a system similar (but not identical, 
as molecular connections may vary in individual pathways) to 
TLR9 network which signals through MyD88 but not through 
TRIF or a condition similar to TLR3 which signals through TRIF 
but not through MyD88 respectively, then in both cases we found 
that ERK = 1 and also both positive and negative paths exist. This 
doesn't allow us to distinguish TLR receptor specific nature of 
wiring between ERK-1/2 and IL-10. Under such conditions one 
can consider building TLR receptor specific models, conduct high 
throughput experiments, train the models with the high through- 
put data [9] and then compute the signaling paths. 

Finally, we tested the model's ability to reproduce experimental 
results and demonstrated its potential therapeutic utility. Several 
knockout conditions were implemented in the model and 
simulations considering such perturbations corroborated with the 
experiments. Similarly, we inhibited ( = 0) the network elements 
that correspond to parasitic infections where our simulations 
reproduced the experimentally reported inhibition of cytokines. 
The therapeutic utility of the model was examined by predicting 
the potential drug targets through MIS analysis. In addition, 
SHP1 -mediated regulation of TLR4 signaling was simulated by 
adding three additional reactions in our model as SHP1 inhibits 
TLR4-mediated IRAKI activation [47,48]] regulating the pro- 
duction of various cytokines. But we did not find any changes in 



the feedback loop mediated regulation of TLR4 signaling (data not 
shown) as SHP-1 acts in the receptor level [47] and it fulfils the 
definition of an inhibitor but not of a feedback element. 

Taken together, this study helps in understanding the structural 
organization based working principles governing the whole cell 
level TLR4 signaling and poses a number of conjectures that can 
be experimentally verified. As proposed earlier [6,9], the model 
can further be trained with high throughput experimental data; 
that will contribute towards enhancing the predictive power of the 
model. 

Materials and Methods 

For qualitative analysis of the TLR4 signal transduction 
pathway, we have used a Boolean logic based modeling 
approaches for modeling signaling networks, implemented in the 
MATLAB toolbox CellNetAnalyzer [6]. The first step in the 
construction of the model was extraction of the information from 
the comprehensive map of the TLR signaling [8], followed by 
incorporation of newer information from literature [Table SI]. 
The information extracted from the comprehensive map of the 
TLR was converted into the "Logical Interaction Hypregraph" 
(LIH) [6]. In the logical interaction hypergraph the interaction 
among the species are the combination of the three simple 
Boolean functions AND, OR and NOT gates or their case specific 
linear combinations. Figure 2 shows the conversion of stoichio- 
metric knowledge to 4 different types of Boolean relation [AND, 
OR, NOT gates] , with a representative set of reactions from the 
total of 263 reactions of the network. 
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Table 1. MISs for inhibiting production of the IL-10 in Yersinia infection. 



S. No. Species 

1 IKK_a.b.g 

2 !c-Fos 

3 ITPL2 

4 NF-k-b_p65.50 

5 TTP 

6 IERK1.2 

7 PP-2A.B 

8 PP2A 

9 HL-10 

10 MKP3 

11 !MKK1 

12 !ERK1.2_Nu 

13 !c-Fos_P 

14 NF-kB.p50 IlkBa IlkBb 

15 ESCIT NF-kB.p50 !lkBb 

16 NF-kB.p50 IlkBb TRAF3 

17 NF-kB.p50 IlkBb TBK1 

18 MEKK1 NF-kB.p50 IlkBb 

19 ILIGAND Ubcl 3.Uev1 A.TIFA.TRAF6.IRAK1 NF-kB.p50 IlkBb 

20 ILIGAND Ubc13.Uev1A.TIFA.TRAF6.IRAK1* NF-kB.p50 IlkBb 

21 Ubc13.Uev1A.TIFA.TRAF6.IRAK1 IMKK3 NF-kB.p50 IlkBb 

22 Ubc13.Uev1A.TIFA.TRAF6.IRAK1* IMKK3 NF-kB.p50 IlkBb 

23 Ubc13.UevlA.TIFA.TRAF6.IRAK1 IMAPKAPK2 NF-kB.p50 IlkBb 

24 Ubc13.Uev1A.TIFA.TRAF6.IRAK1* IMAPKAPK2 NF-kB.p50 IlkBb 

! Denote the permanent inactivation of the species and the species without I sign are species which should be permanently activated for fulfilling target goal. 
* Shows activated form of the IRAKI complex. 
doi:1 0.1 371 /journal.pone.0092481 .t001 



For example, TLR4 receptor and ligand binding steps are 
represented by two interaction graphs (Figure 2; 2 nd row, 2 nd 
column). Biologically it is possible to form a ligand receptor 
complex [TLR4.L] only when both ligand and receptor comes 
together. Using the Boolean gate AND (Figure 2; 2 nd row, 3 rd 
column) we could implement such biological relations where the 
function AND allows the TLR4.L to be active (state = 1) only 
when TLR4 = 1 and LIGAND = 1 . 

On the contrary, there are situations where any one activator (of 
several activators) is sufficient for triggering activation of a 
substrate. For example Mnkl can be phosphorylated either by 
ERK-1/2 or p38MAPK. So (Figure 2; 3 rd row, 2 nd column) the 
OR gate is appropriate to capture the situation (table 1; 3 rd row, 
3 rd column). When any one of ERK-1/2 or p38MAPK= 1, the 
activation state of Mnkl = 1. 

In another scenario, elements like Elkl requires p38MAPK= 1 
and PP2B = 0 for its activation state to become 1 (Figure 2; 5 
row, 1 st column). Boolean representation of such situations is 
shown (Figure 2; 5 th row, 3 rd column) where AND NOT gate 
captured the biological relation between these three species. 

After the conversion of such molecular interactions into the 
logical hypergraphs, the network was constructed for its pic- 
torial representation using the software CellDesigner [49]. The 
picture was exported from CellDesigner to CellNetAnalyzer for 
analysis. 



Computation of the feedback loops in phosphatases 
knockout condition 

The signaling pathways activated through the TLR4 involved 
the activation of the various phosphatases of the system. Role of 
these phosphatases in regulation of feedback loops is calculated in 
the 16 different knockout situations. To knockout a phosphatase 
from the system, it was excluded during the feedback calculations 
by assigning it a fixed value = 0. 

Computation of the dependency matrix 

The dependency matrix of the full-scale network is extremely 
large (181x181) making the analysis process complicated. 
However the information of influence of one species on rest of 
the model elements could be extracted from the dependency 
matrix using CNA [6]. Thus, we selected a set of kinases; 
representatives of various pathways activated during the TLR4 
signaling and analyzed their interdependencies in the knockout 
conditions of the phosphatases. 

Computation of minimal intervention set (MIS) 

Minimal intervention set (MIS) in an interaction networks is the 
(minimal) sets of elements that are to be removed or to be added in 
order to achieve a certain goal such as activation or inhibition of a 
target molecule [6] . MISs were calculated in infection condition to 
abrogate IL-10 production by assigning 0 value to IL- 1 0 node and 
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declaring the fixed activation state for TLR4, Ligand and other 
nodes that are required for receptor-lignad complex formation. To 
minimize the computation time, flag was set for the nodes that 
were not to be included during MIS calculation, like receptor and 
ligand. 

Logical Steady State Analysis 

Logical steady state shows the activation state of a set of 
molecules on assigned time scale for given input signal. For 
computation of the logical steady state in a particular situation, for 
example, in a Myd88 knockout situation, values = 1 or 0 were 
assigned to the respective species and then the LSS were computed 
at time scale 7.1. Time scale 7.1 was selected for analyzing the 
effect of a knockout or addition of inhibitors in the system because 
at time scale 7.1 the signal was allowed to flow through the entire 
span of the network. During the calculation of the LSS few 
interactions were not considered. For example, p38MAPK and 
TAK.TAB1 which are activated at the same time scale are 
connected though a NOT gate, blocking the downstream flow of 
signal. As the p38aMAPK activation and the TAK.TAB1 
activation occur in the time scale 2, and as p38aMAPK inhibits 
the TAK.TAB 1 , calculation of LSS at any further time scale above 
2 calculates the downstream signaling effect considering TAK.- 
TAB 1 = 0. Hence the LSS calculations were carried out without 
including the feedback loop from p38MAPK to TAK.TAB 1. A 
total of five interactions (out of total 263 interactions) were 
removed during the calculation of logical steady state(shown in 
Table S9). 

Requirement of two models for interaction graph and 
hypergraphs based analysis 

As discussed above in the results sections, activation of 
inhibitory loops on the same time scale would result in the 
blockage of the signal through the network, thus the inhibitory 
loops must be removed to allow the signal flow through the 
network. However, removal of the interactions also results in the 
change in the topological properties of the network such as the 
number of feedback loops. To overcome these problem two 
models were built, one model for the logical steady state analysis 
where five reactions were not considered (mentioned above). The 
other model was used for analysis of topological properties such as 
feedback loops or the interdependencies among the species which 
contains all the interactions. Nevertheless, simulated logical steady 
states devoid of the five reactions corroborated with the 
experimental reports suggesting that the signal flow in the model 
is preserved in general and it can be used for making biologically 
testable predictions. 

Model validation using experimentally reported 
scenarios 

We implemented results of various experiments in the logical 
model aiming to validate its predictive power. For example, 
experimental data shows that in Myd88 knockout mice production 
of inflammatory cytokine is down regulated. We knocked out 
My88 in the model by assigning it a 0 value and computed the 
LSS at time scale 7.1 which qualitatively reproduced the 
experimental observations. Several such examples are shown. 
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